#this file estimates the effect at mock thresholds
rm(list=ls(all=TRUE))
graphics.off()
library(rdd)
library(xtable)

setwd() # set working directory
load("data/data_alt_boot.rdata")
#drop all from closed and split lists 
data_alt <- data_alt[which(data_alt$marg.los!=0 & data_alt$marg.win != 1 & 
                             data_alt$openlist == 1 & abs(data_alt$margsim) < 1
                           & data_alt$splitlist == 0 & 
                             data_alt$simscore > 0 & data_alt$simscore < 1),]
#if cases are extremely close they might have wrong signed simscore. Recode these to arbitrarily small simscore
data_alt$margsim[data_alt$margsim <= 0 & data_alt$elected == 1] <-  0.001
data_alt$margsim[data_alt$margsim >= 0 & data_alt$elected == 0] <- -0.001
data_alt <- na.omit(data.frame(data_alt[c("margsim", "electedt1", "rerun")]))

## Placebotests. 

placest <- matrix(NA, nrow=11, ncol=2)
placse  <- matrix(NA, nrow=11, ncol=2)
loop_vec <- seq(from = -5, to = 5, by = 1)/10
for (j in seq(1:11)){
  i <- loop_vec[j]
  placmod       <- RDestimate(electedt1 ~ margsim, data = data_alt, cutpoint = i)
  placmod1      <- RDestimate(rerun     ~ margsim, data = data_alt, cutpoint = i)
  placest[j,1]  <- placmod$est[1]
  placest[j,2]  <- placmod1$est[1]
  placse[j,1]   <- placmod$se[1]
  placse[j,2]   <- placmod1$se[1]
}

#create plot
graphics.off()

quartz(w=8,h = 4)
par(mar=c(4,4,2,1), mfrow = c(1,2))
adj <- 2.5
plot(placest[,1], loop_vec*-10, 
     xlim = c(-.5,.5),
     ylim = c(-5,5),
     pch = 19,
     axes = FALSE,
     ylab = "Cutpoint",
     xlab = "Estimate of rerunning and elected at cutpoint")
axis(1, 
     at = c(-0.5,0,0.5))

par(las = 1)

axis(2,
     at = loop_vec*-10,
     tick = FALSE,
     labels = round(loop_vec,1))


segments(placest[,1] + qnorm(0.025)*placse[,1], 
         loop_vec*-10,
         placest[,1] + qnorm(0.975)*placse[,1], 
         loop_vec*-10,
         lwd = 2)

lines(c(0,0),c(-7.5,7.5),
      lty = 2,
      col = "grey50")

plot(placest[,2], loop_vec*-10, 
     xlim = c(-.5,.5),
     ylim = c(-5,5),
     pch = 19,
     axes = FALSE,
     ylab = "",
     xlab = "Estimate of rerunning at cutpoint")
axis(1, 
     at = c(-0.5,0,0.5))
segments(placest[,2] + qnorm(0.025)*placse[,2], 
         loop_vec*-10,
         placest[,2] + qnorm(0.975)*placse[,2], 
         loop_vec*-10,
         lwd = 2)


lines(c(0,0),c(-7.5,7.5),
      lty = 2,
      col = "grey50")


## Put estimates in table for SI
si_table <- matrix(NA, nrow = nrow(placest), ncol = 4)

for (i in 1:nrow(placest)){
  si_table[i,c(1,3)] <- round(c(placest[i,]),3)
  si_table[i,c(2,4)] <- paste("(", round(c(placse[i,]) ,3), ")", sep="")
}

colnames(si_table) <- c("p(rerunning and elected)"," ", "p(rerunning", "  ")
rownames(si_table) <- loop_vec

xtable(si_table)
